% Simple Model BGRP 
% Jordan Roulleau-Pasdeloup
% Lausanne, September 12, 2015

% This code computes the Ramsey optimal fiscal/monetary policy
% in the baseline New Keynesian model with public capital

%----------------------------------------------------------------
% 1. Declaration
%----------------------------------------------------------------


var c n mc pi phi1 phi2 phi3 phi4 xi i gc exppi rrate;

varexo shock_xi shock_g;

parameters g, phi_r, gamma, sigmaa,  kappa, beta,
delta, phi_pi, b, alpha, rho_xi, Nstar, Kstar, Ucstar, Unstar, Uccstar, Gcstar, Ugcstar, Uggcstar,
Ucnstar, Unnstar, psi, theta, Cstar, T, chi_g, Gstar, Ystar, Ynstar, Ykstar, Ynnstar, Ynkstar, Ykkstar, rho_g;

%----------------------------------------------------------------
% 2. Parameters
%----------------------------------------------------------------

g = 0.2;
theta = 6;
gamma = 0.2857;
sigmaa = 2;
kappa = 0.03;
psi = 90;

T=16;
beta = 0.99;
b=0.0;
delta = 0.02;
phi_pi = 1.5;
rho_xi  = 0.9;
rho_g = .9;

alpha = b*beta^T/(1-beta*(1-delta))*delta/g;

Nstar = gamma/( (1-g)*(1-gamma)+gamma );

Cstar = ((1-g)*Nstar)^(1/(1-b))*(alpha*g/(delta*(1-g)))^(b/(1-b));

Kstar = alpha/delta*g/(1-g)*Cstar;

Ystar = Nstar*Kstar^b;

Ynstar = Kstar^b;

Ykstar = b*Nstar*Kstar^(b-1);

Ynnstar = 0;

Ynkstar = b*Kstar^(b-1);

Ykkstar = b*(b-1)*Nstar*Kstar^(b-2);

Gstar= g*Ystar;

Gcstar = Gstar - delta*Kstar;

Ucstar = gamma*Cstar^(gamma*(1-sigmaa)-1)*(1-Nstar)^((1-gamma)*(1-sigmaa));

Unstar = -(1-gamma)*Cstar^(gamma*(1-sigmaa))*(1-Nstar)^((1-gamma)*(1-sigmaa)-1);

Uccstar = -(1+gamma*(sigmaa-1))*Ucstar/Cstar;

Ucnstar = (1-gamma)*(sigmaa-1)*Ucstar/(1-Nstar);

Unnstar = -( (1+(1-gamma)*(sigmaa-1))*(1-gamma)*Cstar^(gamma*(1-sigmaa))*(1-Nstar)^((1-gamma)*(1-sigmaa)-2) );

Ugcstar = Gcstar^-sigmaa;

Uggcstar = -sigmaa*Gcstar^(-sigmaa-1);

phi2star = 0;

phi3star = Ucstar;

chi_g = phi3star*((1-alpha)*Gstar)^sigmaa;

%----------------------------------------------------------------
% 3. Model
%----------------------------------------------------------------

model (linear);

Nstar*Unnstar/Unstar*n + Cstar*Ucnstar/Unstar*c + xi =  phi3 + (Ynnstar*Nstar/Ynstar)*n + theta/(Ynstar*Unstar)*((Unnstar*Ucstar-Ucnstar*Unstar)/Ucstar^2)*phi2
                                                        -theta*(Unstar/Ucstar)*(Ynnstar*Nstar/Ynstar^2)*phi2
                                                        -phi1(-1)*Ucnstar/(beta*Unstar) + Ucnstar*phi1/Unstar; // FOC wrt N
                          
Cstar*Uccstar/Ucstar*c + Nstar*Ucnstar/Ucstar*n +xi =  theta/(Ynstar*Ucstar)*((Ucnstar*Ucstar-Uccstar*Unstar)/Ucstar^2)*phi2 + phi3
                                                        -phi1(-1)*Uccstar/(Ucstar*beta) + phi1*Uccstar/Ucstar; // FOC wrt C

pi = beta*pi(+1) + (theta-1)/psi*mc; // New Keynesian Phillips Curve

mc = (-Uccstar*Cstar/Ucstar + Ucnstar*Cstar/Unstar)*c + (Nstar*Unnstar/Unstar-Nstar*Ucnstar/Ucstar-Ynnstar*Nstar/Ynstar)*n; // Definition of marginal cost

c = c(+1)+Ucstar/(Cstar*Uccstar)*(i-pi(+1)+log(beta)-(1-rho_xi)*xi)+Nstar*Ucnstar/(Cstar*Uccstar)*(n(+1)-n); // Euler Equation

Ynstar*Nstar/Ystar*n = (1-g)*c + (1-alpha)*g*gc;

exppi = pi(+1);

xi = rho_xi*xi(-1)-shock_xi;

phi4 = -beta*Ucstar*phi1;

Ystar*Ucstar*pi + phi2 = + phi2(-1)-0.000000001*phi1(-1);

phi1 = 0.000000001*phi4;

gc = rho_g*gc(-1)+shock_g;

rrate = i-pi(+1)+log(beta);

end;

%----------------------------------------------------------------
% 4. Computation
%----------------------------------------------------------------

initval;
i=-log(beta);
end;
steady;
resid(1);

check;



shocks; 
var shock_xi;
stderr 0.1;
var shock_g;
stderr 0.01;
end; 


stoch_simul(irf=0,periods=1000);

break;

figure(3);
subplot(331);
plot((1-g)*c_shock_xi/abs(xi_shock_xi(1)),'-','LineWidth',3,'color','blue');
grid on
title('Private Consumption')
subplot(332);
plot(alpha*g*gi_shock_xi/abs(xi_shock_xi(1)),'-','LineWidth',3,'color','blue');
grid on
title('Government Investment')
subplot(333);
plot(n_shock_xi/abs(xi_shock_xi(1)),'-','LineWidth',3,'color','blue');
grid on
title('Hours Worked')
subplot(334);
plot(b*beta^T/(1-beta*(1-delta))*k_shock_xi/abs(xi_shock_xi(1)),'-','LineWidth',3,'color','blue');
grid on
title('Public Capital')
subplot(335);
plot(xi_shock_xi/abs(xi_shock_xi(1)),'-','LineWidth',3,'color','blue');
grid on
title('Preference Shock')
subplot(336);
plot(i_shock_xi+1/beta-1,'-','LineWidth',3,'color','blue');
grid on
title('Nominal Rate')
subplot(337);
plot(pi_shock_xi/abs(xi_shock_xi(1)),'-','LineWidth',3,'color','blue');
grid on
title('Inflation Rate')
subplot(338);
plot(mc_shock_xi/abs(xi_shock_xi(1)),'-','LineWidth',3,'color','blue');
grid on
title('Marginal Cost')

print(gcf,'-dpsc2','Ramsey_NK_noF');


Util=zeros(100,1);

for i=1:100;
    
    if sigmaa>1
    
    Util(i) = beta^(i-1) * exp(xi_shock_xi(i)) * (((Cstar*(1+c_shock_xi(i)))^gamma*(1-Nstar*(1+n_shock_xi(i)))^(1-gamma))^(1-sigmaa)-1)/(1-sigmaa)
            + beta^(i-1) * exp(xi_shock_xi(i)) * chi_g*(Gstar*(1+gc_shock_xi(i))^(1-sigmaa)-1)/(1-sigmaa);
            
    else
        
    Util(i) = beta^(i-1) * exp(xi_shock_xi(i)) * (gamma*log((Cstar*(1+c_shock_xi(i))))+(1-gamma)*log(1-Nstar*(1+n_shock_xi(i))))        
            + beta^(i-1) * exp(xi_shock_xi(i)) * chi_g*log(Gstar*(1+gc_shock_xi(i)));
    
    end
end;

Welf_NK = sum(real(Util))
